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2 . ABSTRACT 

"^ ■ We use two very large cosmological simulations to study how the density profiles of 

C^ \ relaxed ACDM dark halos depend on redshift and on halo mass. We confirm that 

these profiles deviate slightly but systematically from the NFW form and are better 
approximated by the empirical formula, dlogp/dlogr ex r Q , first used by Einasto to fit 
star counts in the Milky Way. The best-fit value of the additional shape parameter, a, 
increases gradually with mass, from a ~ 0.16 for present-day galaxy halos to a ~ 0.3 
for the rarest and most massive clusters. Halo concentrations depend only weakly 

r-^. , on mass at z = 0, and this dependence weakens further at earlier times. At z ~ 3 

C^> ' the average concentration of relaxed halos does not vary appreciably over the mass 

range accessible to our simulations (M > 3 x 10 11 h~ 1 M & ). Furthermore, in our biggest 
simulation, the average concentration of the most massive, relaxed halos is constant at 
(C200) ~ 3-5 to 4 for < z < 3. These results agree well with those of Zhao et al (2003b) 

f~*) and support the idea that halo densities reflect the density of the universe at the time 

they formed, as proposed by Navarro, Frenk & White (1997). With their original 

parameters, the NFW prescription overpredicts halo concentrations at high redshift. 

This shortcoming can be reduced by modifying the definition of halo formation time, 

jJJ ■ although the evolution of the concentrations of Milky Way mass halos is still not 

reproduced well. In contrast, the much- used revisions of the NFW prescription by 
Bullock et al. (2001) and Eke, Navarro & Steinmetz (2001) predict a steeper drop 
in concentration at the highest masses and stronger evolution with redshift than are 
compatible with our numerical data. Modifying the parameters of these models can 
reduce the discrepancy at high masses, but the overly rapid redshift evolution remains. 
These results have important implications for currently planned surveys of distant 
clusters. 

Key words: methods: N-body simulations - methods: numerical -dark matter - 
galaxies: haloes - galaxies:structure 



1 INTRODUCTION the linear power spectrum from which nonlinear structures 

have grown. As a result, it is useful to parametrize halo pro- 
Over the past decade, cosmoloeical N-body simulations have ci l ■ l ■ • l r l l ,, . ii 

° , tiles by simple empirical formulae, such as that proposed by 

shown consistently that equilibrium dark matter halos have AT „ , „ „,, .. ,.,„„_ .„„„ .,„„„ , c, , Tr ,„ n 

, . „ J , , „, Navarro, Frenk & White (1995, 1996, 1997, hereafter NFW): 

spherically-averaged mass density prohles which are approx- 
imately "universal" in form; i.e., their shape is independent 
of mass, of the values of the cosmological parameters, and of 
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where p cr it = 3H 2 /SivG is the critical density for closurqj, S c 
is a characteristic density contrast, and r s is a scale radius. 
Note that this formula contains two scale parameters but no 
adjustable shape parameter. 

As discussed in some detail by NFW and confirmed by 
subsequent numerical work, the two parameters of the NFW 
profile do not take arbitrary values, but are instead corre- 
lated in a way that reflects the mass-dependence of halo 
assembly times (e.g., Kravtsov, Klypin & Khokhlov 1997; 
Avila-Reese et al. f999; Jing 2000; Ghigna et al. 2000; Bul- 
lock et al. 200f ; Eke, Navarro & Steinmetz 2001; Klypin et 
al. 2001). The basic idea behind this interpretation is that 
the characteristic density of a halo tracks the mean density 
of the universe at the time of its formation. Thus, the later 
a halo is assembled, the lower its characteristic density, <5 C , 
or, equivalently, its "concentration" (see Section 12.21 for a 
definition) . 

Although the general validity of these trends is well 
established, a definitive account of the redshift and mass 
dependence of halo concentration is still lacking, even for 
the current concordance cosmology. This is especially true 
at high masses, where enormous simulation volumes are re- 
quired in order to collect statistically significant samples of 
these rare systems. Simulating large cosmological volumes 
with good mass resolution is a major computational chal- 
lenge, and until recently our understanding of the mass pro- 
file of massive halos has been rather limited, derived largely 
from small numbers of individual realizations or from ex- 
trapolation of models calibrated on different mass scales 
(NFW; Moore et al. 1998; Ghigna et al. 2000; Klypin et al. 
2001; Navarro et al. 2004; Tasitsiomi et al. 2004; Diemand 
et al. 2004; Reed et al. 2005). 

Individual halo simulations may result in biased con- 
centration estimates, depending on the specific selection cri- 
teria used to set them up. In addition, they are unlikely 
to capture the full scatter resulting from the rich variety 
of possible halo formation histories. Extrapolation based on 
poorly tested models can also produce substantial errors, as 
recently demonstrated by Neto et al. (2007). These authors 
analyzed the mass-concentration relation for halos identi- 
fied at z — in the Millennium Simulation of Springel et al. 
(2005, hereafter MS) and confirmed the earlier conclusion of 
Zhao et al. (2003b) that the models of Bullock et al. (2001, 
hereafter B01) and Eke et al. (2001, hereafter ENS) (which 
were calibrated to match galaxy-sized halos) severely un- 
derestimate the average concentration of massive clusters, 
by up to a factor of ~ 3. 

Estimates of concentrations can also be biased by the 
inclusion of unrelaxed halos. These often have irregular den- 
sity profiles caused by major substructures. Smooth density 
profiles are often poor fits to such halos, and the resulting 
concentration estimates are ill-defined because they depend 
on the radial range of the fit and choice of weighting. They 
can also lead to spurious correlations (see for example fig- 
ure 9 of Neto et al. 2007). Consequently, in this paper we 
follow Neto et al. and select only relaxed halos for analy- 
sis. This is not without its own problems. Such selection 
biases against recently formed halos, which may preferen- 



1 We express Hubble's constant as H(z) and its present-day value 
as H(z = 0) = H = 100 h km s _1 Mpc" 1 . 



tially have lower concentrations. However we believe this is 
preferable to polluting the sample with meaningless concen- 
tration estimates of the kind that arise when smooth spheri- 
cal models are force-fit to lumpy, multi-modal mass distribu- 
tions. Hayashi & White (2008) stacked all halos, regardless 
of dynamical state, in the MS and studied the resulting mean 
profiles as a function of halo mass. The relatively small dif- 
ferences between their results and those found below shows 
that the inclusion of unrelaxed halos has rather little effect 
in the mean. 

A further preoccupation concerns indications that halo 
profiles deviate slightly but systematically from the NFW 
model (Moore et al. 1998, Jing & Suto 2000, Fukushige & 
Makino 2001, Navarro et al. 2004, Prada et al. 2006, Merritt 
et al. 2006), raising the possibility that estimating concen- 
trations by force-fitting simple formulae to numerical data 
may result in subtle biases that could mask the real trends. 
This is especially important because of hints that such de- 
viations depend systematically on halo mass (Navarro et al. 
2004, Merritt et al. 2005) . Evaluating and correcting for such 
deviations is important in order to establish conclusively the 
mass and redshift dependence of halo concentration. 

These uncertainties are unfortunate since observations, 
especially at high redshift, often focus on exceptional sys- 
tems. For example, massive galaxy clusters are readily iden- 
tified in large-scale surveys of the distant universe, and un- 
derstanding their internal structure will be critical for the 
correct interpretation of cluster surveys intending to con- 
strain the nature of dark energy. These will make precise 
measurements of the evolution of cluster abundance in sam- 
ples detected by gravitational lensing, by their optical or X- 
ray emission, or through the Sunyaev-Zel'dovich effect (see, 
e.g., Carlstrom et al. 2002, Hu 2003, Majumdar & Mohr 
2003, Holder 2006, and references therein). 

There is at present no ab initio theory that can reliably 
predict the internal structure of CDM halos. The models of 
Zhao et al (2003a), Wechsler et al (2002) (see also Lu et 
al. 2006) are interesting, but as shown in Neto et al (2007) 
they account for only a small fraction in the measured dis- 
persion in concentration at a fixed mass. There is currently 
no substitute for direct numerical simulation when detailed 
predictions are needed for comparison with observation. 

We address these issues here by combining results from 
the MS with results from an additional simulation which fol- 
lowed a substantially smaller volume but with better mass 
resolution. This allows us to extend the range of halo masses 
for which we can measure concentrations and to assess how 
these measures are affected by numerical resolution. Our 
analysis procedure follows closely that of Neto et al. (2007) . 
In particular, we concentrate in this paper on the properties 
of halos which are relaxed according to the criteria defined 
by these authors; mean density profiles for all MS halos of 
given mass, regardless of dynamical state, are presented by 
Hayashi & White (2008). We begin in Section by describ- 
ing briefly the numerical simulations and the halo catalogue 
on which this study is based. In Section [3] we present our 
main results for the dependence of profile shape and con- 
centration on halo mass and redshift. We conclude with a 
brief discussion and summary in Section [4] 
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Figure 1. Left panels: The stacked density profile of 464 halos of virial mass ~ 2x 10 14 h _1 MQ, identified at z = in the MS. The 
curves in the upper left panel show different NFW fits obtained by varying the radial range fitted: [0.02, l]r v i r (solid red), [0.05, l]r v i r 
(dashed green) and [0.1,l]r v j r (dot-dashed blue). Note the small disagreement between the actual profile shape and the NFW model 
(upper panels). This leads to concentration estimates which depend slightly on the innermost radius of the fit. Fits using eq. ((4} are 
more robust to such variations in fitting range, as shown in the lower left panel. Panels on the right show the residuals corresponding to 
the fits shown on the left. 



2 THE NUMERICAL SIMULATIONS 

The analysis presented here is based primarily on halos iden- 
tified in the Millennium Simulation (MS) of Springel et al. 
(2005). The halo identification and cataloguing procedure 
follows closely that described in detail by Neto et al. (2007). 
For completeness, we here recapitulate the main aspects of 
the procedure, referring the interested reader to the earlier 
papers for details. 



2.1 Simulations 

The MS is a large N-body simulation of the concordance 
ACDM cosmogony. It follows N = 2160 particles in a peri- 
odic box of side Lbox = 500 h~ 1 Mpc. The chosen cosmolog- 
ical parameters were fi m = Sldm + fib — 0.25, fib = 0.045, 
h = 0.73, fi A = 0.75, n = 1, and a$ = 0.9. Here fi de- 
notes the present-day contribution of each component to the 
matter-energy density of the Universe, expressed in units of 
the critical density for closure, p C rit; n is the spectral in- 
dex of the primordial density fluctuations, and as is the rms 



linear mass fluctuation in a sphere of radius 8h 1 Mpc at 
2 = 0. The particle mass in the MS is 8.6 x 1O 8 /i _1 M . 
Particle interactions are softened on scales smaller than the 
(Plummer-equivalent) softening length, e = 5/i -1 kpc. 

We also use a second simulation of a smaller volume 
(100 3 h~ 3 Mpc 3 ) within the same cosmological model. This 
simulation followed iV = 900 3 particles of mass 9.5 x 
10 7 /i _1 Mq and softened interactions on scales smaller than 
e = 2.4/i _1 kpc. It thus has approximately 9 times better 
mass resolution than the MS. Hereafter we refer to it as the 
hMS. 



2.2 Halo Catalogues 

Our halo cataloguing procedure starts from a standard 
b — 0.2 friends-of-friends list of particle groups (Davis et 
al. 1985) and refines it with the help of SUBFIND, the sub- 
halo finder algorithm described by Springel et al. (2001). 
Each halo is centred at the location of the minimum of the 
gravitational potential of its main SUBFIND subhalo, and 
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Figure 2. Left panel: The best-fit shape parameter, a (eq. J4)l), as a function of halo mass and redshift, after binning and stacking halos 
by mass. Solid and dashed lines correspond to the MS and hMS simulations, respectively. Only results corresponding to halos with more 
than ~ 500 particles and to stacks with more than 10 halos are shown. Note the good general agreement between the two simulations. 
Differences are substantial only where the number of particles is less than 3000. Right panel: Values of a for halos with more than 
3000 particles plotted as a function of the dimensionless "peak height" parameter u(M,z), defined as the ratio between the critical 
overdensity <5 cr it(z) for collapse at redshift z and the linear rms fluctuation at z within spheres containing mass M. The larger the value 
of v the rarer and more massive the corresponding halo. Note that this scaling accounts satisfactorily for the redshift dependence of the 
mass-concentration relation shown in the left-hand panel. 



this centre is used to compute the virial radius and mass of 
each halo. 

We define the virial radius, ta, as that of a sphere of 
mean density A x p cx w- Note that this defines implicitly the 
"virial mass" of the halo, Ma, as that enclosed within ta- 
The choice of A varies in the literature. The most popular 
choices are: (i) a fixed value, as in NFW's original work, 
where A = 200 was adopted; or (ii) a value motivated by 
the spherical top-hat collapse model, where A ~ 178 f2 n ; 45 
for a fiat universe (Eke et al. 1996). The latter choice gives 
A = 95.4 at 2 = for the cosmological parameters adopted 
for the MS. We keep track of both definitions in our halo 
catalogue, and will specify our choice by a subscript. Thus, 
M200 and T2oo are the mass and radius of a halo adopting 
A = 200, whereas M v i r and r vlI are the values corresponding 
to A = 95.4 at z — 0. Quantities specified without subscript 
assume A = 200, so that, e.g., M — Mwa. The halo con- 
centration is defined as the ratio between the virial radius 
and the scale radius: c = C200 = r2oo/r 3 . In this case, the 
characteristic density is related to the concentration by 



5 C = (200/3) c 3 /[ln(l + c) - c/(l + c)]. 



(2) 



Since halos are dynamically evolving structures, we use 
a combination of diagnostics in order to flag non-equilibrium 
systems. Following Neto et al. (2007), we assess the equilib- 
rium state of each halo by measuring: (i) the substruc- 
ture mass fraction; i.e., the total mass fraction in re- 
solved substructures whose centres lie within r2oo, /sub = 
y"l,-Jn b Msub.i /A/200; (ii) the centre of mass displace- 
ment, s = |r c — r cm |/r2oo, defined as the normalized offset 
between the location of the minimum of the potential and 
the barycentre of the mass within r2oo; and (iii) the virial 
ratio, 2T/\U\, where T is the total kinetic energy of the halo 
particles within r v j r and U their gravitational self-potential 
energy. 



Using these criteria, we shall consider halos to be re- 
laxed if they satisfy all the following conditions: / su b < 0.1, 
s < 0.07, and 2T/\U\ < 1.35. (See Neto et al. 2007 for full 
details.) We shall also impose a minimum number of parti- 
cles in order to be able to say something meaningful about 
internal halo structure. We initially set this limit at 500 par- 
ticles, but following the analysis of profile shapes presented 
in Section 3.2 we subsequently adopt a stricter criterion of 
3000 particles. We consider only relaxed halos in this study, 
since only for such systems can the mass profiles of indi- 
vidual objects be represented accurately by a simple fitting 
formula with a small number of parameters. Such formulae 
are also useful for characterizing the average profiles of large 
ensembles of halos, since the fluctuations in the individual 
systems then average out. Hayashi & White (2007) present 
such mean profiles as a function of mass for all MS halos re- 
gardless of their dynamical state. We compare their results 
with our own below. 

With these restrictions, including the 3000 particle 
limit, our final MS catalogue contains 128233, 77190, 30603, 
and 9392 relaxed halos at z — 0, 1, 2, and 3, respectively. 
(The corresponding numbers for the hMS catalogue are 8131, 
6652, 4112 and 2194.) The overall fractions of these halos 
that are relaxed in the MS are 78%, 60%, 50%, and 48% 
at these redshifts, respectively. We note that these frac- 
tions also depend on halo mass: at z — about ~ 85% 
of ~ 10 12 /i _1 M0 halos are relaxed by our criteria, but only 
~ 60% of ~ 10 15 /i _1 M Q halos. In order to obtain usefully 
large statistical samples, we restrict our analysis to the red- 
shift range < z < 3 in the following. 



3 HALO DENSITY PROFILES 

For each halo in the samples described in Section [2] we have 
computed a spherically-averaged density profile by measur- 
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ing the halo mass in 32 equal intervals of log 10 (r) over the 
range > log 10 (r/r vir ) > -2.5. 

Profiles may also be stacked in order to obtain an aver- 
age profile for halos of similar mass. This procedure has the 
advantage of erasing individual deviations from a smooth 
profile which are typically due to the presence of substruc- 
ture. Such deviations increase the (already considerable) 
scatter in the parameters fitted to individual profiles and 
may mask underlying trends in the data. The left panels in 
Fig. [l] show the profile that results from stacking 464 ha- 
los of mass ~ 2 x 10 14 /i _1 Mq identified at z — in the 
MS. We choose to show r 2 p vs r rather than p vs r in order 
to remove the main radial trend and enhance the dynamic 
range of the plot. Similar stacked profiles for all halos in the 
MS (regardless of dynamical state) are shown by Hayashi & 
White (2008). 



3.1 Deviations from NFW and concentration 
estimates 

The concentration of a halo is defined using the scale radius 
of the profile, r 3 . This identifies the location of the maxi- 
mum of the r 2 p profile. However, as shown in Fig. \T\ the 
peak is rather broad, leading to some uncertainty in its ex- 
act location when noise is present. Typically this problem 
is addressed by fitting the numerical data to some specified 
functional form over an extended radial range. The curves 
in the top-left panel of Fig. [l] show the result of fitting the 
stacked halo profile with the NFW formula (eq. fTJ), but 
varying the radial range of the fit as shown by the labels. 
This results in slightly different estimates for r s , and conse- 
quently for the concentration, c v i r = Jvir/fV Increasing the 
innermost radius of the fit from 0.02 to O.lrVir results in a 
concentration estimate that decreases from ~ 7.5 to ~ 6.7. 

This variation is a result of the slight (but significant) 
mismatch between the shape of the NFW profile and that of 
the stacked halo, as shown by the residuals in the top-right 
panel of Fig. [T] The "S" shape of the residuals implies that 
the stacked profile steepens more gradually with radius in 
a log-log plot than does the NFW profile, a result that has 
been discussed in detail by Navarro et al. (2004), Prada et 
al. (2006) and Merritt et al. (2006). 

These results suggest that force-fitting NFW profiles 
may induce spurious correlations between mass and concen- 
tration. In particular, when halos are identified in a sin- 
gle cosmological simulation, the numerical resolution varies 
systematically with halo mass (less massive halos are more 
poorly resolved) so that the radial range available for fitting 
is a strong function of halo mass. 

One way to address this issue is to adopt a fitting for- 
mula that better matches the mean profile of simulated ha- 
los. As discussed by Navarro et al. (2004), Prada et al. (2006) 
and Merritt et al. (2006), improved fits are obtained adopt- 
ing a radial density law where the logarithmic slope is a 
power-law of radius, 



dlogp 
dlogr 



= -2 



(3) 




which implies a density law of the form 

In(p/p- 3 ) = -(2/a)[(r/r- 3 )°'-l]. (4) 

Here p_2 is the density at r = r_2. This density law was 
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Figure 3. Three stacked halo density profiles with different values 
of v, and at different redshifts, as labelled. The profiles are chosen 
to illustrate the variation in halo structure as a function of v 
indicated in Fig. 2. The larger the halo mass, the larger the value 
of a and the sharper the curvature in the density profile as a 
function of radius. 



first introduced by Einasto (1965) who used it to describe 
the distribution of old stars within the Milky Way. For con- 
venience, we will thus refer to it as the Einasto profile. Note 
that according to our definition, r_2 corresponds to the ra- 
dius where the logarithmic slope of the density profile has 
the "isothermal" value, —2. In this sense, r_2 is equivalent 
to the NFW scale-length, r s , and again marks the location 
of the maximum in the r 2 p profile shown in Fig. [T] 

The bottom- left panel of Fig. [T] shows that, at the cost 
of introducing an extra shape parameter, the fits improve 
to the point that the sensitivity of concentration estimates 
to the fitted radial range is effectively eliminated. Thus, 
concentrations obtained by fitting eq. @ to the stacked 
halo profiles are robust against variations in fitting details. 
Hayashi & White (2007) show that the same is true for fits 
to stacks of all halos of a given mass, rather than just the 
relaxed halos used to make Fig. [TJ and indeed, the a and c 
values they find are very similar to the values we obtain here, 
showing that our restriction to relaxed halos has relatively 
little effect in the mean. 

We conclude that accounting for the subtle difference 
between halo profile shape and the NFW fitting formula is 
worthwhile in order to avoid possible fitting-induced biases 
in concentration estimates. In the remainder of this paper 
we quote concentrations, ca = »"a/?~_2, which are estimated 
by fitting density profiles by eq. @. We discuss in the next 
section how a is chosen for these fits. 



3.2 The mass and redshift dependence of profile 
shape 

The above discussion suggests that the shape parameter, 
a, should be used to improve the description of the typical 
density profiles of simulated halos and to eliminate possible 
biases in estimates of their concentration. To do this, it is 
necessary to understand whether (and how) a varies with 
redshift and/or halo mass. 
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Figure 4. Distributions of concentration parameters estimated 
using eq. Q for the 464 individual halo profiles stacked in Fig. [T] 
The solid (black) histogram corresponds to fits where a was ad- 
justed separately for each individual halo, the dashed (red) his- 
togram to fits where a was set to the value implied by the a(v) 
relation of Fig. [2] (eq. J5}). The numbers in the legend give the 
median and scatter of the distributions. The excellent agreement 
between the two distributions indicates that unbiased and accu- 
rate concentration estimates for relaxed halos may be obtained 
by fixing a to the value determined by eq. l5l l. 



We explore this in Fig. [2] The left panel shows how 
the best-fit value of a depends on halo mass for the average 
profiles of relaxed halos stacked according to their mass. We 
consider only halos with at least 500 particles within the 
virial radius, and stacks containing at least fO halos. The 
solid and dashed curves in this plot correspond to the MS 
and hMS simulations, respectively. 

This figure illustrates several interesting points. In the 
first place, we note that there is good agreement between the 
two simulations for halos which are represented by at least 
3000 particles, but that systematic differences are visible 
when the MS halos contain fewer particles than this. In the 
rest of this paper we will thus only present results for halos 
containing at least 3000 particles within the virial radius. 

A second interesting point is that there are well-defined 
trends for a to increase with mass at each redshift, and 
with redshift at each mass. The (weak) trend with mass was 
already visible in Navarro et al. (2004) and Merritt et al. 
(2005), although the small number of halos in these studies 
made their results rather inconclusive in the face of the large 
halo-to-halo scatter. The use of stacked profiles in Fig. [2] 
together with the much larger number of halos in the simu- 
lations used here, leads to a far more convincing demonstra- 
tion of the trends than was previously possible. 

The trend with redshift at given mass may seem sur- 
prising, but its interpretation is made clear by the right 
panel of Fig. [2] Here we show a as a function of a di- 
mensionless "peak-height" parameter, defined as v(M, z) = 
8 CT it(z)/a(M, z); i.e., as the ratio of the linear density thresh- 
old for collapse at z to the rms linear density fluctuation 
at z within spheres of mean enclosed mass M. The pa- 
rameter v(M, z) is related to the abundance of objects of 
mass M at redshift z. v(M*,z) — 1 defines the charac- 
teristic mass, M*(z), of the halo mass distribution at red- 



shift z. v(M, z) 3> 1 then corresponds to rare objects with 
M 3> M* , while v(M, z) <C 1 corresponds to objects in the 
low-mass tail of the distribution. The parameter v plays a 
key role in the Press-Schechter model for the growth of non- 
linear structure (see, for example, Lacey & Cole 1993). 

The right panel of Fig. [5] shows that all curves coincide 
when expressed in terms of v(M, z) (and when considering 
halos with N > 3000 particles). Thus, the redshift depen- 
dence in the left panel merely reflects the fact that objects 
of given virial mass correspond to very different values of v 
at different redshifts. It is interesting that the dependence 
of a on v is very weak for v < 1 (it is nearly constant at 
a ~ 0.16), but it increases rapidly for rarer, more massive 
objects, reaching a ~ 0.3 for v ~ 3.5. A simple formula, 



a = 0.155 + 0.0095 u 



(5) 



describes the numerical results quite accurately. 

Taylor & Navarro (2001), Austin et al (2005) and 
Dehnen & McLaughlin (2005) investigated a halo model 
based on the assumption that the phase space density, p/er' 5 , 
was a simple powerlaw of radius. Interestingly, the density 
profile they found is almost identical to an Einasto profile 
with 0.14 < a < 0.18 and so the sharp upturn we find in a 
at v > 1 could be taken as indication that such a model is 
not valid for rare and massive halos. 

The dependence of profile shape on v is illustrated in 
Fig. [3] where we plot r 2 p profiles for halo stacks correspond- 
ing to three different values of v. 1.0, 2.0, and 3.2. The larger 
v, the larger a, and the more sharply the profile peaks. It 
is unclear what causes these systematic trends, but the fact 
that they depend on v (rather than, say, on mass alone) is 
an important clue for models that attempt to explain the 
near-universality of dark halo density profiles. While our re- 
sults here are based purely on relaxed halos, very similar 
results were found by Hayashi & White (2007) for the aver- 
age profiles of stacks of all MS halos, regardless of dynamical 
state. 



3.3 Concentration estimates 

As we discussed above, a fitting formula other than NFW is 
needed to obtain concentration estimates that are insensitive 
to the radial range fitted. Einasto's profile, eq. Q, accom- 
plishes this at the cost of introducing an additional shape pa- 
rameter, a. Adjusting a to fit the detailed structure of indi- 
vidual halos would negate the spirit of the NFW programme 
which attempts to predict the structure of dark halos in any 
hierarchical cosmology from its initial power spectrum and 
the global cosmological parameters. Three-parameter fits to 
individual halos are also susceptible to strongly correlated 
parameter errors. We have therefore explored whether fixing 
q as a function of halo mass and redshift according to the 
a{v) relation of eq. <(Sj results in significantly different con- 
centration distributions than adjusting it freely to fit each 
halo. 

We show the result in Fig. |4j which shows concentra- 
tion distributions for the 464 individual halos which were 
stacked to make Fig. [l] We compare the distribution ob- 
tained from full three-parameter Einasto fits to that ob- 
tained when q is set to the value predicted for each halo 
from its mass. We find that for most halos the fits give es- 
sentially the same concentration, and the distributions of 
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Figure 5. The mass dependence of concentration as a function of redshift. Concentration estimates are derived from Einasto fits to the 
density profile of "relaxed" halos over the radial range [0.05, l]r v i r . Points, boxes and whiskers show the median and the 5, 25, 75, and 
95 percentiles of the concentration distribution within each mass bin. Black symbols correspond to the MS, red symbols to the smaller, 
but higher resolution hMS. Only results for halos with more than 3000 particles are shown for each simulation. The solid-dotted, dashed, 
dot-dashed and solid lines show the models of NFW, ENS, B01 and the modified NFW, respectively, for comparison with our results. 
The lower dotted line in the upper left panel indicates the relation obtained by Neto et al. (2007) from NFW fits to relaxed halos at 
z = 0. The remaining dotted lines show the powerlaw fits given in Table 1. 



concentration are almost indistinguishable. The medians of 
the fixed-a and floating-a distributions are 6.85 and 6.90, re- 
spectively. The scatters also coincide, as noted in the labels 
of Fig. 13] Interestingly, the scatter around the best-fit profile 
is typically only slightly smaller for the Einasto model than 
for the NFW model, showing that the difference between the 
two models is much smaller than the deviation of the profile 
of typical individual halos from either model. In this sense, 
the 3-parameter Einasto profile has little advantage over the 
2-parameter NFW profile when estimating dark matter halo 
concentrations. 

This is easily understood. Estimating the concentration 
of a halo is equivalent to locating the "peak" of the r 2 p 
profile. Provided that the shape of the fitted profile is, on 
average, a good approximation to the simulated ones, no sys- 
tematic offset is expected between concentrations estimated 
using the two procedures. We conclude that concentrations 
may be estimated robustly by fitting Einasto profiles to indi- 
vidual halos with a fixed to the value obtained from eq. ©. 
All values reported below were obtained using this proce- 
dure. 



3.4 The mass and redshift dependence of halo 
concentration 

The concentration-mass relation is shown in Fig. [S] for 
z = 0, 1, 2 and 3. The upper-left panel is equivalent to Fig. 6 
of Neto et al. (2007), except that our concentrations are es- 
timated from fits covering the range [0.05, l]r v ir using an 
Einasto rather than an NFW profile, and we only show re- 
sults for halos with more than 3000 particles. The differences 
are small, as may be judged from the power-law fit proposed 
by Neto et al, c 20 o = 5.26 (M2oo/lO 14 ft" 1 M )" 010 , shown 
as a dotted line in the upper-left panel of Fig. [5] This power 
law fit is very similar (after correction for differing concen- 
tration and mass definitions) to that which Maccio et al. 
(2007) obtained from NFW fits to halos in an independent 
set of (smaller) simulations, despite the fact that these au- 
thors included halos with as few as 250 particles, which we 
would consider to be significantly under-resolved on the ba- 
sis of our own tests. 

Both the concentration values and the trends with mass 
and redshift are very similar to those presented in figure 2 
of Zhao et al (2003b) who analysed a set of ACDM Simula- 
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tions of varying resolution. The small offsets between their 
mean concentrations and our results are consistent with the 
slight differences we have noted when switching from NFW 
to Einasto models for determining concentrations. In ad- 
dition to confirming these earlier results, the much larger 
volume of our simulations results in a better determination 
of the intrinsic scatter about the mean relation. 

Red and black symbols in Fig. [5] correspond to the MS 
and hMS simulations, respectively, with dots plotted at the 
median concentration for each mass bin. Boxes represent the 
lower and upper quartiles of the concentration distributions, 
while the whiskers show their 5% and 95% tails. We also 
show the concentrations predicted by three previously pro- 
posed analytic models: NFW (solid-dotted magenta), B01 
(dot-dashed black), ENS (dashed green). As discussed by 
Neto et al. (2007), none of these models reproduces the sim- 
ulation results over the full mass range accessible at z — 0: 
the NFW and ENS models underestimate the concentration 
of galaxy-sized halos, whereas the B01 model fails dramati- 
cally at the high-mass end, where it predicts a sharp decline 
which is not seen in the simulations, [j 

There is a hint in the z — panel that the relation 
is flattening at the high-mass end, where the NFW predic- 
tions at z = appear slightly better than those of ENS. This 
is because a constant concentration for very rare and mas- 
sive objects is implicit in the NFW model, which assumes 
that the characteristic density of a halo reflects that of the 
universe at the time it collapsed. Very massive systems as- 
sembled very recently, and therefore share the same collapse 
time (i.e., they are being assembled today). 

The near-constant concentration of the most massive 
halos is considerably more obvious at higher redshift, pre- 
sumably because well-resolved halos (i.e., those with N > 
3000 in the MS or hMS) become rarer and rarer with increas- 
ing lookback time. Indeed, whereas at z = our MS halo 
catalogue spans the range 0.75 < v < 3, at z = 3 all the 
halos retained have v>3. As a result, the average concen- 
tration at z — 3 is almost independent of mass over the 
accessible mass range, i.e., M>3 x lO 11 ft _1 M0. It is inter- 
esting that the average concentration of the most massive 
halos (i.e. v > 3) is similar at all redshifts, c ~ 3.5 to 4. 
This evokes the proposals of Zhao et al. (2003a,b), Tasit- 
siomi et al. (2004) and Lu et al. (2006), all of whom argue 
that halos undergoing rapid growth should all have similar 
concentration. 

The evolution of the mass-concentration relation seen 
in our numerical simulations is not predicted by any of the 
published prescriptions. The original NFW model shows a 



2 We note that NFW, ENS and B01 parametrize the initial 
power spectrum in slightly different ways, and that the pre- 
dicted concentrations are sensitive to the exact choice of pa- 
rameters. For example, NFW and ENS use the parameter T to 
characterize the shape of the linear power spectrum. We use 
r = 0.15 here since that provides the best fit to the actual 
power spectrum used in the MS. For B01 we have used the de- 
fault values advocated in the latest version of their software 
(K = 2.8 and F = 0.001 in their notation), which is available 
from http://www.physics.uci.edu/~bullock/CVIR/ The differ- 
ences between the predictions shown here and in Fig. 4 of Hayashi 
& White (2008), or in Neto et al (2007), for example, are due to 
slight differences in the values chosen for these parameters. 



Table 1. Values of the constants A and B for the best straight- 
line fit J6j to the data shown in Fig. \5\ 



Redshift 


A 


B 





-0.138 


2.646 


0.5 


-0.125 


2.372 


1 


-0.092 


1.891 


2 


-0.031 


0.985 


3 


-0.004 


0.577 



flattening of the concentration-mass relation with increas- 
ing redshift, similar to that observed in the simulations, but 
it predicts insufficient evolution in concentration at given 
mass. As a result, this model overestimates the concen- 
trations by an increasing amount with increasing redshift, 
about 40% at z = 3. The ENS and B01 models fail to 
reproduce the features of the mass-concentration-redshift 
relations at high mass, predicting a stronger mass depen- 
dence and much more evolution than is seen. In these two 
models, the concentration of halos of given mass scales in- 
versely as the expansion factor, so that shape of the mass- 
concentration relation remains fixed. While the high-mass 
discrepancy between the B01 model and our measurements 
can be reduced by changing the parameters to K = 2.8 
and F = 0.001 (see Wechsler et al 2006), neither for this 
model nor for the ENS model can parameter changes pro- 
duce agreement with the weak redshift evolution seen at 
high mass both here and by Zhao et al. (2003b). On the 
other hand, the ENS and B01 models predict the concen- 
tration evolution of galaxy mass halos substantially better 
than the NFW model. Note that our simulation data do 
not disagree significantly from those analysed by ENS and 
B01. The discrepancies result from extrapolation of their 
proposed relations beyond the range where they were reli- 
ably tested. 

In the NFW model, the definition of formation time in- 
volves two physical parameters, F and /, and a proportion- 
ality constant, C, that relates the value of the characteristic 
halo density to the mean density of the universe at the time 
of collapse. The formation redshift of a halo is taken to be 
the redshift at which a fraction, F, of its mass was first con- 
tained in progenitors each individually containing at least 
a fraction / of its mass. In the original NFW prescription, 
F = 0.5, / = 0.01 and C = 3000. We find that the ob- 
served trend in the slope of the mass-concentration relation 
with redshift can be approximately reproduced by simply 
changing F to F — 0.1, keeping / = 0.01 as before; the nor- 
malization of the relation is then approximately reproduced 
by taking C = 600. The resulting curves are shown as solid 
blue lines in Fig. [5] This modified NFW model succeeds well 
in matching the redshift evolution, but its mass dependence 
is still too shallow at z — and too steep at z = 3, leaving 
room for improvement, especially at low masses where the 
B01 and ENS prescriptions give better predictions of the 
evolution rate. We have checked that the same model also 
gives an acceptable fit for other cosmologies, in particular for 
a simulation of a ACDM model similar to hMS but with the 
values of the cosmological parameters inferred from the 3- 
year WMAP satellite data (Spergel et al. 2007) (fi = 0.236, 
fi A = 0.764, n — 0.97 and a 8 — 0.74). We have also checked 
that this modified NFW model gives a good description of 
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the concentrations found in the scale-free models presented 
in Navarro, Frenk & White (1997) which span a range of 
spectral indices and have either £1 — 1 or £1 — 0.1. If more 
accurate results for the particular cosmology assumed in this 
paper are desired, we provide the coefficients of power-law 
fits of the form 



log 10 (c 2 oo) = Alog 10 (M 2 oo) + B 



(6) 



in Table [T] and plot these best fit curves in Fig. [5] 

The disagreement between the simulation data and the 
original NFW prescription is up to a few tens of percent. 
The BOl and ENS prescriptions underestimate the concen- 
tration of the most massive halos by factors of ~ 2 to ~ 3 
at high redshift. Since the characteristic density of a halo 
scales roughly as the third power of the concentration (see 
eq. ©), this implies that the characteristic densities of such 
halos are underestimated by at least an order of magnitude 
by the latter two models, leading to dramatic changes in 
the expected lensing power, X-ray luminosity, and S-Z de- 
tectability of massive clusters at high redshift. 

It is also interesting to compare the z = 
concentration-mass relation found here to the one which 
Hayashi & White (2007) obtained by fitting Einasto profiles 
to stacks of all halos of a given mass, rather than to stacks 
of relaxed halos. The results in their Fig. 4 lie very close 
to the NFW relation plotted in the upper-left panel of our 
own Fig. [S] Thus, the restriction to relaxed halos has little 
systematic effect on Einasto-based concentrations at masses 
above about 3 x 10 13 /i _1 Mq, but at lower masses it results 
in somewhat larger concentrations. This differs slightly from 
the result of Neto et al (2007) who found z = concentra- 
tions obtained by fitting NFW (rather than Einasto) pro- 
files to relaxed halos to be greater by a larger amount at 
high mass. The median of their NFW-based concentrations 
is also closer to the NFW model prediction on galaxy scales 
(both for relaxed and for all halos) than is the median of the 
Einasto-based concentrations which we plot in Fig. [5] This 
is at least in part a result of the bias illustrated in Fig. Q] 

Neto et al. (2007), Hayashi & White (2007) and this pa- 
per are all based on the same simulations (although Hayashi 
& White do not use the MS). The small but (statistically) 
significant differences in their derived concentration-mass re- 
lations show that at the 10 to 20% level these relations de- 
pend on the details of the halo sample selection and profile 
fitting procedures. The modified NFW model works rea- 
sonably well over the full range of mass and redshift stud- 
ied here, certainly rather better than the revisions proposed 
by ENS and BOl. Significant discrepancies remain, however, 
notably at galaxy masses where the NFW model underpre- 
dicts the rate of evolution, resulting in systematically low 
Einasto-based concentrations at z ~ and systematically 
high values at z > 2. For such objects the evolution rate 
is better predicted by the ENS and BOl models. Given the 
good agreement between various simulations on this mass 
scale (see, e.g., Zhao et al 2003b, Maccio et al. 2007 and 
Neto et al. 2007), we believe these concentration estimates 
for low redshift ACDM galaxy halos to be accurate. This 
implies that the NFW model needs further improvement for 
such objects. At higher masses and higher redshifts, how- 
ever, the (modified) NFW formalism works quite well, and 
should be preferred to those of BOl or ENS. 

With our modified parameters, the NFW prescription 



not only fits concentrations in the WMAP1 cosmology in- 
vestigated here, but also our (less extensive) set of data for 
the WMAP3 cosmology. Without detailed testing, however, 
it is unclear if the prescription can be extended to other 
regimes of interest. For example, what concentrations are 
expected for dwarf galaxy halos or for "first object" halos at 
z ~ 30? The incorrect conclusions reached by applying the 
BOl and ENS formulae to massive clusters at high redshift 
are testimony to the dangers of using empirical formulae 
outside the range where they have been tested. Clearly, fur- 
ther theoretical effort aimed at understanding the factors 
which determine the internal structure of dark halos would 
be a welcome complement to the numerical work presented 
here. 



4 SUMMARY 

We have used data from the Millennium Simulation and 
from a smaller but higher resolution simulation to examine 
how the density profiles of relaxed ACDM halos vary with 
halo mass and redshift. We study profile shape and con- 
centration for halos with mass exceeding ~3x lQ ll h~ Mq 
over the redshift range < z < 3. Our main results may be 
summarized as follows. 

• We confirm the conclusion of previous studies that, in 
the mean, the shape of spherically averaged ACDM halo 
density profiles deviates slightly but systematically from the 
two-parameter fitting formula proposed by Navarro, Frenk 
& White (1995, 1996, 1997). A more accurate description 
is provided by the three-parameter Einasto (1965) profile, 
for which the logarithmic slope is a power-law of radius, 
dlogp/dlogr oc r a . We show that this fitting formula gives 
concentration estimates which are insensitive to the radial 
range fitted, albeit at the price of an additional shape pa- 
rameter. Although Einasto fits avoid small but significant 
biases that arise if all halos are fit to the NFW model, we 
emphasise that individual halos typically deviate from either 
model by more than the difference between them. 

• Using stacked profiles of halos of similar mass, we show 
that the shape parameter, a, depends systematically on halo 
mass and on redshift. These dependences can be collapsed 
into a dependence on the single "peak height" parameter, 
v{M,z) = 5 CI it{z)/a(M,z). Halos with large v are rare ob- 
jects in the high-mass tail of the halo mass distribution and 
have large a values (see eq. (©). This provides an important 
clue for models that attempt to interpret the dependence of 
halo density profiles on mass and redshift. 

• The dependence of halo concentration on mass becomes 
progressively weaker with increasing redshift, as found ear- 
lier by Zhao et al (2003b). At Z — 3, concentrations are 
almost independent of mass over the mass range accessi- 
ble to our simulations. Relaxed halos with v > 3 (the 
rarest and most massive systems) have similar concentra- 
tions, (C200) = 3.5 to 4, at all the redshifts we have studied 

• The models of Bullock et al. (2001) and Eke et al. (2001) 
fail to reproduce our measured concentrations for high-mass 
and high-redshift objects, predicting a stronger mass depen- 
dence and more evolution than is seen in the simulations. 
Parameters in these models can be changed to reduce the 
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strength of their mass dependence, but the predicted red- 
shift evolution, while fitting galaxy mass halos well up to 
redshift z — 1, remains substantially too strong at high mass 
and at higher redshifts. As a result, the predictions of these 
models for high-redshift galaxy clusters can be in error by 
large factors. The original model of Navarro et al. (1997) 
overpredicts the concentrations of such objects at redshifts 
beyond 1 (by up to ~ 40% at z = 3) but a modified NFW 
model with a different definition of formation redshift re- 
produces the simulation results substantially better over the 
redshift and mass ranges we have examined. Both the orig- 
inal and the modified NFW models underestimate the con- 
centration evolution of relaxed 10 12 /i _1 Mq halos, leading to 
30% discrepancies at z = and z = 3. 

We hope that our simulation results will stimulate the- 
oretical work aimed at a deeper understanding of the factors 
which determine the internal structure of ACDM halos. Such 
work may eventually result in simple recipes like those of 
NFW, ENS and B01. Substantial errors are found, however, 
when these published prescriptions are extrapolated beyond 
the regime where their authors tested them, in particular, 
to the regime relevant to high-redshift galaxy clusters. This 
demonstrates that careful numerical work is mandatory be- 
fore any recipe can be applied in a new regime. When making 
forecasts for surveys of distant massive clusters, our simula- 
tions show our modified NFW recipe to give reliable results 
at least out to Z — 3. 
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